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Abstract 

We develop a moment-equation-copula-closure method for the inexpensive approximation 
of the steady state statistical structure of strongly nonlinear systems which are subjected to 
correlated excitations. Our approach relies on the derivation of moment equations that describe 
the dynamics governing the two-time statistics. These are combined with a non-Gaussian pdf 
representation for the joint response-excitation statistics, based on copula functions that has 
i) single time statistical structure consistent with the analytical solutions of the Fokker-Planck 
equation, and ii) two-time statistical structure with Gaussian characteristics. Through the 
adopted pdf representation, we derive a closure scheme which we formulate in terms of a con¬ 
sistency condition involving the second order statistics of the response, the closure constraint. 
A similar condition, the dynamics constraint^ is also derived directly through the moment equa¬ 
tions. These two constraints are formulated as a low-dimensional minimization problem with 
respect to the unknown parameters of the representation, the minimization of which imposes 
an interplay between the dynamics and the adopted closure. The new method allows for the 
semi-analytical representation of the two-time, non-Gaussian structure of the solution as well 
as the joint statistical structure of the response-excitation over different time instants. We 
demonstrate its effectiveness through the application on bistable nonlinear single-degree-of- 
freedom energy harvesters with mechanical and electromagnetic damping, and we show that 
the results compare favorably with direct Monte-Garlo simulations. 


1 Introduction 

In numerous systems in engineering, uncertainty in the dynamics is as important as the known 
conservation laws. Such an uncertainty can be introduced by external stochastic excitations, e.g. 
energy harvesters or structural systems subjected to ocean waves, wind excitations, earthquakes, 
and impact loads [T]-[^. For these cases, deterministic models cannot capture or even describe 
the essential features of the response and to this end, understanding of the system dynamics and 
optimization of its parameters for the desired performance is a challenging task. On the other 
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hand, a probabilistic perspective can, in principle, provide such information but then the challenge 
is the numerical treatment of the resulted descriptive equations, which are normally associated with 
prohibitive computational cost. 

The focal point of this work is the development of a semi-analytical method for the inexpen¬ 
sive probabilistic description of nonlinear vibrational systems of low to moderate dimensionality 
subjected to correlated inputs. Depending on the system dimensionality and its dynamical char¬ 
acteristics, numerous techniques have been developed to quantify the response statistics, i.e. the 
probability density function (pdf) for the system state. For systems subjected to white noise, 
Fokker-Planck-Kolmogorov (FPK) equation provides a complete statistical description of the re¬ 
sponse statistics [7]-[^. However, exact analytical solutions of the FPK equation are available only 
for a small class of systems. An alternative computational approach, the path integral solution 
(PIS) method, has been developed to provide the response pdf for general nonlinear systems at a 
specific time instant given the pdf of an earlier time instant. Many studies have been focused on the 
application of step-by-step PIS method numerically 10 -12 and analytically 13 15 reporting its 


effectiveness on capturing the response statistics. On the other hand, for non-Markovian systems 
subjected to correlated excitations the joint response-excitation pdf method provides a computa¬ 
tional framework for the full statistical solution 16-18 . However, such methodologies rely on the 


solution of transport equations for the pdf and they are associated with very high computational 
cost especially when it comes to the optimization of system parameters. 

To avoid solving the transport equations for the pdf, semi-analytical approximative approaches 
with significantly reduced computational cost have been developed. Among them the most popular 
method in the context of structural systems is the statistical linearization method H^ 23 , which 


can also handle correlated excitations. The basic concept of this approach is to replace the original 
nonlinear equation of motion with a linear equation, which can be treated analytically, by minimiz¬ 
ing the statistical difference between those two equations. Statistical linearization performs very 
well for systems with unimodal statistics, i.e. close to Gaussian. However, when the response is 
essentially nonlinear, e.g. as it is the case for a double-well oscillator, the application of statistical 
linearization is less straightforward and involves the ad-hoc selection of shape parameters for the 
response statistics 
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An alternative class of methods relies on the derivation of moment equations, which describes 
the evolution of the the joint response-excitation statistical moments or (depending on the nature 
of the stochastic excitation) the response statistical moments [^ 27 . The challenge with moment 
equations arises if the equation of motion of the system contains nonlinear terms in which case 
we have the well known closure problem. This requires the adoption of closure schemes, which 
essentially truncate the infinite system of moment equations to a finite one. The simplest approach 
along this line is the Gaussian closure 28 but nonlinear closure schemes have also been developed 
(see e.g. 29-^). In most cases, these nonlinear approaches may offer some improvement compared 
with the stochastic linearization approach applied to nonlinear systems but the associated compu¬ 
tational cost is considerably larger 38 . For strongly nonlinear systems, such as bistable systems, 
these improvements can be very small. Bistable systems, whose potential functions have bimodal 
shapes, have become very popular in energy harvesting applications [39}|4^, where there is a need 


for fast and reliable calculations that will be able to resolve the underlying nonlinear dynamics in 
order to provide with optimal parameters of operation (see e.g. 


47 M). 


The goal of this work is the development of a closure methodology that can overcome the 
limitations of traditional closure schemes and can approximate the steady state statistical structure 
of bistable systems excited by correlated noise. We first formulate the moment equations for the 


2 






















joint pdf of the response and the excitation at two arbitrary time instants 49 . To close the 


resulted system of moment equations, we formulate a two-time representation of the joint response- 
excitation pdf using copula functions. We choose the representation so that the single time statistics 
are consistent in form with the Fokker-Planck-Kolmogorov solution in steady state, while the joint 
statistical structure between two different time instants is represented with a Gaussian copula 
density. Based on these two ingredients (dynamical information expressed as moment equations 
and assumed form of the response statistics), we formulate a minimization problem with respect 
to the unknown parameters of the pdf representation so that both the moment equations and the 
closure induced by the representation are optimally satisfied. For the case of unimodal systems, 
the described approach reproduces the statistical linearization method while for bi-modal systems 
it still provides meaningful and accurate results with very low computational cost. 

The developed approach allows for the inexpensive and accurate approximation of the second 
order statistics of the system even for oscillators associated with double-well potentials. In addition, 
it allows for the semi-analytical approximation of the full non-Gaussian joint response-excitation 
pdf in a post-processing manner. We illustrate the developed approach through nonlinear single- 
degree-of-freedom energy harvesters with double-well potentials subjected to correlated noise with 
Pier son-Moskowitz power spectral density. We also consider the case of bi-stable oscillators coupled 
with electromechanical energy harvesters (one and a half degrees-of-freedom systems), and we 
demonstrate how the proposed probabilistic framework can be used for performance optimization 
and parameters selection. 


2 Description of the Method 

In this section, we give a detailed description of the proposed method for the inexpensive compu¬ 
tation of the response statistics for dynamical systems subjected to colored noise excitation. The 
computational approach relies on two basic ingredients: 

• Two-time statistical moment equations. These equations will be derived directly from the 
system equation and they will express the dynamics that govern the two-time statistics. For 
systems excited by white-noise, single time statistics are sufficient to describe the response 
but for correlated excitation, this is not the case and it is essential to consider higher order 
moments. Note that higher (than two) order statistical moment equations may be used but 
in the context of this work two-time statistics would be sufficient. 

• Probability density function (pdf) representation for the joint response-excitation statistics. 
This will be a family of probability density functions with embedded statistical properties 
such as multi-modality, tail decay properties, correlation structure between response and 
excitation, or others. The joint statistical structure will be represented using copula functions. 
We will use representations inspired by the analytical solutions of the dynamical system when 
this is excited by white noise. These representations will reflect features of the Hamiltonian 
structure of the system and will be used to derive appropriate closure schemes that will be 
combined with the moment equations. 

Based on these two ingredients, we will formulate a minimization problem with respect to the 
unknown parameters of the pdf representation so that both the moment equations and the closure 
induced by the representation are optimally satisfied. We will see that for the case of unimodal 
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systems the described approach reproduces the statistical linearization method while for bi-modal 
systems it still provides meaningful and accurate results with very low computational cost. 

For the sake of simplicity, we will present our method through a specific system involving a 
nonlinear SDOF oscillator with a double well potential. This system has been studied extensively 
in the context of energy harvesting especially for the case of white noise excitation 41 50 52 


However, for realistic setups it is important to be able to optimize/predict its statistical properties 
under general (colored) excitation. More specifically we consider a nonlinear harvester of the form 


X \x kix -h ksx^ = y. 


( 1 ) 


where x is the relative displacement between the harvester mass and the base, y is the base excita¬ 
tion representing a stationary stochastic process, A is normalized (with respect to mass) damping 
coefficient, and ki and ks are normalized stiffness coefficients. 


y(t) 


h(t) 



Figure 1: Nonlinear energy harvester with normalized system parameters. 


2.1 Two-time Moment System 

We consider two generic time instants, t and s. The two-time moment equations have been consid¬ 
ered previously in for the determination of the solution of a ‘half’ degree-of-freedom nonlinear 
oscillator by utilizing a Gaussian closure. We multiply the equation of motion at time t with the 
response displacement x{s) and apply the mean value operator □ (ensemble average). This will 
give us an equation which contains an unknown term on the right hand side. To determine this 
term we repeat the same step but we multiply the equation of motion with y{s). This gives us the 
following two-time moment equations: 

x{t)y{s) + \x(t)y{s) + kix{t)y{s) + k 3 x{t)^y{s) = y{t)y{s), (2) 

x{t)x{s) + Xx{t)x{s) + kix{t)x{s) + ksx{t)^x{s) = y{t)x{s). (3) 


Here the excitation is assumed to be a stationary stochastic process with zero mean and a given 
power spectral density; this can have an arbitrary form, e.g. monochromatic, colored, or white 
noise. Since the system is characterized by an odd restoring force, we expect that its response also 
has zero mean. Moreover, we assume that after an initial transient the system will be reaching 
a statistical steady state given the stationary character of the excitation. Based on properties of 
mean square calculus , we interchange the differentiation and the mean value operators. Then 
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the moment equations will take the form: 


dt^' 

- 


\-l 

dt“ 

d- 



(4) 

^y{t)x{s). 

(5) 


Expressing everything in terms of the covariance functions, above equations will result in: 


d 


d 


d‘^ 


+ \-C% + k,C% + k,x{tfy{s) = 

+ x^cli + k,cli + ksWMs) = ^cli, 

where the covariance function is defined as 

Cil = x(t)y{s) = C^yit -s)= C^yij). 


(6) 

(7) 

(8) 


Taking into account the assumption for a stationary response (after the system has gone through 
an initial transient phase), the above moment equations can be rewritten in terms of the time 
difference r = t — s: 


d‘^ 


^-^Cxx{ 


02 

(9) 


(10) 


Note that all the linear terms in the original equation of motion are expressed in terms of covariance 
functions, while the nonlinear (cubic) terms show up in the form of fourth order moments. To 
compute the latter we will need to adopt an appropriate closure scheme. 


2.2 Two-time PDF representations and induced closures 

In the absence of higher-than-two order moments, the response statistics can be analytically ob¬ 
tained in a straightforward manner. However, for higher order terms it is necessary to adopt an 
appropriate closure scheme that closes the infinite system of moment equations. A standard ap¬ 
proach in this case, which performs very well for unimodal systems, is the application of Gaussian 


closure which utilizes Isserlis’ Theorem 53 to connect the higher order moments with the second 


order statistical quantities. Despite its success for unimodal systems, Gaussian closure does not 
provide accurate results for bistable systems. This is because in this case (i.e. bistable oscilla¬ 
tors) the closure induced by the Gaussian assumption does not reflect the properties of the system 
attractor in the statistical steady state. 

Here we aim to solve this problem by proposing a non-Gaussian representation for the joint 
response-response pdf at two different time instants and for the joint response-excitation pdf at two 
different time instants. These representations will: 

• incorporate specific properties or information about the response pdf (single time statistics) 
in the statistical steady state. 
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• capture the correlation structure between the statistics of the response and/or excitation at 
different time instants by employing Gaussian copula density functions, 


• have a consistent marginal with the excitation pdf (for the case of the joint response-excitation 
pdf). 


2.2.1 Representation Properties for Single Time Statistics 


We begin by introducing the pdf properties for the single time statistics. The selected representation 
will be based on the analytical solutions of the Fokker-Planck equation which are available for the 
case of white noise excitation IE] , and for vibrational systems that has an underlying Hamiltonian 
structure. Here we will leave the energy level of the system as a free parameter - this will be 
determined later. In particular, we will consider the following family of pdf solutions (Figure]^): 


= yex.Y>{-]^U{x)} = ^ exp ■ 




( 11 ) 


where U is the potential energy of the oscillator, 7 is a free parameter connected with the energy 
level of the system, and T is the normalization constant expressed as follows: 



Figure 2 : (a) Representation of the steady state pdf for single time statistics of a system with 
double-well potential. The pdf is shown for different energy levels of the system, (b) The joint 
response excitation pdf is also shown for different values of the correlation parameter c ranging 
from small values (corresponding to large values of |r|) to larger ones (associated with smaller 
values of |r|). 
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2.2.2 Correlation Structure between Two-time Statistics 


Representing the single time statistics is not sufficient since for non-Markovian systems (i.e. corre¬ 
lated excitation) the system dynamics can be effectively expressed only through (at least) two-time 
statistics. To represent the correlation between two different time instants we introduce Gaussian 
copula densities . A copula is a multivariate probability distribution with uniform marginals. 

It has emerged as an useful tool for modeling stochastic dependencies allowing the separation of 
dependence modeling from the given marginals 57 . Based on this formulation we obtain pdf 


representations for the joint response-response and response-excitation at different time instants. 

Joiut respouse-excitatiou pdf. We first formulate the joint response-excitation pdf at two 
different (arbitrary) time instants. In order to design the joint pdf based on the given marginals 
of response and excitation, we utilize a bivariate Gaussian copula whose density can be written as 
follows 


C (u, v) — 


1 




: exp 


2 c$-^ (u) (v) - ($-^ {uf -h (^)^) 


( 13 ) 


where u and v indicate cumulative distribution functions and the standard cumulative distribution 
function is given as the following form: 


^{x) = 



(14) 


Denoting with x the argument that corresponds to the response at time t, with y the argument for 
the excitation at time s = t — r, and with g{y) the (zero-mean) marginal pdf for the excitation, we 
have the expression for the joint response-excitation pdf. 


y) = f{x)g{y)C (F(x), G{y)), 


= f{x)g{y) 


1 


Vl — 


: exp 


2c$-i {F{x)) {G{y)) - ($“' {F{x)f + {G{y)f) 

2(1 -c2) 


( 15 ) 


where c defines the correlation between the response and the excitation and has values — 1 < c < 
1 and F{x) and G{y) are the cumulative distribution functions obtained through the response 
marginal pdf, /(x), and the excitation marginal pdf, g{y), respectively. Note that the coefficient 
c depends on the time difference r = t — s of the response and excitation. This dependence will 
be recovered through the resolved second-order moments (over time) between the response and 
excitation. 

Joint response-response pdf. The joint pdf for two different time instants of the response, 
denoted as p{x,z), is a special case of what has been presented. In order to avoid confusion, a 
different notation ^ is used to represent the response at a different time instant s = t — r. We have: 


p{x, z) = f(x)f(z)C (F(x),F(z )), 


= f{x)f{z) 


Vi -' 


; exp 


2c<I>“^ {F{x)) <I>“^ (F(z)) - ($-1 {F(x)f + (^’( 2 ))^) 

2(1 - c2) 


( 16 ) 


where c is a correlation constant (that depends on the time-difference t). Note that the response 
2 at the second time instant follows the same non-Gaussian pdf corresponding to the single time 
statistics of the response. 
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In Figure [^, we present the above joint pdf with the marginal / (response) having a 
bimodal structure and the marginal g (excitation) having a Gaussian structure. For c = 0 we have 
independence, which essentially expresses the case of very distant two-time statistics, while as we 
increase c the correlation between the two variables increases referring to the case of small values 
of r. 

2.2.3 Induced Non-Gaussian Closures 



Figure 3: The relation between x{t)^x{s) and x{t)x{s). Exact relation is illustrated in red curve 
and approximated relation using non-Gaussian pdf representations is depicted in black curve. 


Using these non-Gaussian pdf representations, we will approximate the fourth order moment terms 
that show up in the moment equations. We numerically observe that in the context of the pdf 
representations given above, the relation between x{t)^x{s) and x{t)x{s) is essentially linear (see 
Figure]^. To this end, we choose a closure of the following form for both the response-response 
and the response-excitation terms: 


x{tfx{s) = x{t)x{s), 


(17) 


where px^x is the closure coefficient for the joint response-response statistics. The value of px^x is 
obtained by expanding both x{t)^x{s) and x{t)x{s) with respect to c keeping up to the first order 
terms: 


'= JJ xzp{x, z)dxdz = 2^ Jxf{x)eTi ^ {2F{x) — 1) dx^ c-hG(c^), (18) 

JJ x^zp{x, z)dxdz = 2^ j x^/(x) erf“^ (2F(x) — 1) dx|| J z/(z) erf“^ (2F(z) — 1) dz|c -h G(c^), 


( 19 ) 


where the error function is given by: 


erf(x) = 



( 20 ) 
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Thus, we observe that the assumed copula function in combination with the marginal densities 
prescribe an explicit dependence between fourth- and second-order moments, expressed through 
the coefficient: 


f x^f(x) erf ^ {2F{x) — 1) dx 
f xf{x) erf“^ {2F{x) — 1) dx 


We emphasize that this elosure eoeffieient does not depend on the time-differenee r but only on the 
single time statisties and in partieular the energy level of the system, defined by 7 . To this end, for 
any given marginal pdf /, we can analytically find what would be the closure coefficient under the 
assumptions of the adopted copula function. 

The corresponding coefficient for the joint response-excitation statistics px,y can be similarly 
obtained through a first order expansion of the moments: 


Px,y 


/ x^f{x) erf ^ {2F{x) — 1) dx 
f xf{x) erf“^ {2F(x) — 1) dx 


( 22 ) 


The closure coefficient px^y has exactly the same form with the closure coefficient px^x and it does 
not depend on the statistical properties of the excitation nor on the time-difference r but only on 


the energy level 7 . We will refer to equations ( 21 ) and ( 22 ) as the elosure eonstraints. This will 


be one of the two sets of constraints that we will include in the minimization procedure for the 
determination of the solution. 


2.2.4 Closed Moment Equations 


The next step involves the application of above closure scheme on the derived two-time moment 
equations. By directly applying the induced closure schemes on equations © and ( [Tq| ), we have 
the linear set of moment equations for the second-order statistics: 


5^ d d‘^ 

'o^^xyi^) + ^-^Cxyix) + {ki + px,yk^)Cxy{T) = (^) 5 

d 

0^2 ^xx F ^~q^Cxx{F^ F {k\ T Px,xk3)Cxx{'^) ~ ^xy ( '^)- 


(23) 

(24) 


Using the Wiener-Khinchin theorem, we transform the above equations to the corresponding power 
spectral density equations: 


{{jujf F Xijuj) FkiF px,yk3}Sxy{^) = {jujfSyy{uj), 
{(j^) k I Px^xks}Sxx{^) — (j^) ^xy{^)- 


(25) 

(26) 


These equations allow us to obtain an expression for the power spectral density of the response 
displacement in terms of the excitation spectrum: 


Sxx{^) — 


UJ 


{ki F Px,yk3 F j{X(j)}{ki -h Px,xk3 - - j{Xuj)} 

Integration of the above equation will give us the variance of the response: 


Syy{uj). 


(27) 


POO no 

— I Sxx{p^')dLo — I 
Jo Jo 


{kl + Px,yk3 - j{Xu)}{ki + Px^xks - j{Xuj)} 


Syy{u)du. (28) 
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The last equation is the second constraint, the dynamics constraint^ which expresses the second order 
dynamics of the system. Our goal is to optimally satisfy it together with the closure constraints 
defined by equations ( 21 ) and ( 22 ). 


2.2.5 Moment Equation Copula Closure (MECC) Method 

The last step is the minimization of the two set of constraints, the closure constraints and the 
dynamics constraint^ which have been expressed in terms of the system response variance x‘^. The 
minimization will be done in terms of the unknown energy level 7 and the closure coefficients px^x 
and Px,y- 

More specifically, we define the following cost function which incorporates our constraints: 


(7? Px,x, Px,y) 



_ ^ ^yy{^) _ 

{ki + px,yk3 - + j{Xoj)}{ki + px,xk3 - 




2 


/x^/(x)erf ^ (2T(x) — 1) dx 1 f /x^/(x)erf ^ (2T(x) — 1) dx 1 

/ xf{x) erf“^ {2F{x) — 1) dx j ykx,y j xf{x) erf“^ {2F{x) — 1) dx J 

(29) 


Note that in the context of statistical linearization only the first constraint is minimized while 
the closure coefficient is the one that follows exactly from a Gaussian representation for the pdf. In 
this context there is no attempt to incorporate in an equal manner the mismatch in the dynamics 
and the pdf representation. The minimization of this cost function essentially allows mismatch 
for the equation (expressed through the dynamic constraint) but also for the pdf representation 
(expressed through the closure constraints). For linear systems and an adopted Gaussian pdf for 
the response the above cost function vanishes identically. 


3 Applications to Bistable Energy Harvesters under Corre¬ 
lated Excitation 

We apply the presented Moment Equation Copula Closure (MECC) method to two nonlinear vibra¬ 
tional systems. One is a single-degree-of-freedom (SDOE) bistable oscillator with linear damping 
that simulates energy harvesting while the other is a similar bistable oscillator coupled with an elec¬ 
tromechanical energy harvester. Eor both applications, it is assumed that the stationary stochastic 
excitation has a power spectral density given by the Pierson-Moskowitz spectrum, which is typical 
for excitation created by random water waves: 

'S'H = Q ^ exp(-T)^ (30) 

where q controls the intensity of the excitation. 

3.1 SDOF Bistable Oscillator Excited by Colored Noise 

Eor the colored noise excitation that we just described, we apply the MECC method. We consider 
a set of system parameters that correspond to a double well potential. Depending on the intensity 
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of the excitation (which is adjusted by the factor q), the response of the bistable system dives’ in 
three possible regimes. If q is very low, the bistable system is trapped in either of the two wells 
while if q is very high the energy level is above the homoclinic orbit and the system performs cross¬ 
well oscillations. Between these two extreme regimes, the stochastic response exhibits combined 
features and characteristics of both energy levels and it has a highly nonlinear, multi-frequency 
character 


58,59 


Despite these challenges, the presented MECC method can inexpensively provide with a very 
good approximation of the system’s statistical characteristics as it is shown in FigureIn particular 
in Figure we present the response variance as the intensity of the excitation varies for two sets 
of the system parameters. We also compare our results with direct Monte-Carlo simulations and 
with a standard Gaussian closure method [^[^|^. 

For the Monte-Carlo simulations the time series for the excitation has been generated as the sum 
of cosines over a range of frequencies. The amplitudes and the range of frequencies are determined 
through the power spectrum while the phases are assumed to be random variables which follow 
a uniform distribution. In the presented examples, the excitation has power spectral density that 
follows the Pier son-Moskowitz spectrum. Once each ensemble time series for the excitation has 
been computed, the governing ordinary differential equation is solved using a 4th/5th order Runge- 
Kutta method. For each realization the system is integrated for a sufficiently long time interval in 
order to guarantee that the response statistics have converged. For each problem, we generate 100 
realizations in order to compute the second-order statistics. However, for the computation of the 
full joint pdf, a significantly larger number of samples is needed reaching the order of 10^. 




Figure 4: Mean square response displacement with respect to the amplification factor of Pierson- 
Moskowitz spectrum for the bistable system with two different sets of system parameters, (a) A = 1, 
ki = —1, and = 1. (b) A = 0.5, ki = —0.5, and k^ = 1. 


We observe that for very large values of q the computed approximation closely follows the Monte- 
Carlo simulation. On the other hand, the Gaussian closure method systematically underestimates 
the variance of the response. For lower intensities of the excitation, the exact (Monte-Carlo) 
variance presents a non-monotonic behavior with respect to q due to the co-existence of the cross- 
and intra-well oscillations. While the Gaussian closure has very poor performance on capturing 
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this trend, the MECC method can still provide a satisfactory approximation of the dynamics. Note 
that the non-smooth transition observed in the MECC curve is due to the fact that for very low 


values of q the minimization of the cost function (29) does not reach a zero value while this is the 


case for larger values oi q. In other words, in the strongly nonlinear regime neither the dynamics 
constraint nor the closure constraint is satisfied exactly, yet this optimal solution provides with a 
good approximation of the system dynamics. 

After we have obtained the unknown parameters 7 , px^x and px^y by minimizing the cost function 
for each given g, we can then compute the covariance functions and the joint pdf in a post-process 
manner. More specifically, since a known 7 corresponds to a specific px,y (equation ( 22 )) we 
can immediately determine Cxyir) by taking the inverse Eourier transform of Sxy found through 
equation (25). The next step is the numerical integration of the closed moment equation ( [2^ 
utilizing the determined value px^x with initial conditions given by 


<^^^( 0 ) = j {x;j)dx, and Cxx{0) = 0, 


(31) 


where the second condition follows from the symmetry properties of Cxx- Note that we integrate 
equation (24) instead of using the inverse Eourier transform as we did for Cxyir) so that we can 


impose the variance found in the last equation by integrating the resulted density for the determined 
7 . Using the correlation functions Cxx{^) and Cxy{r) we can also determine, for each case, the 
correlation coefficient c of the copula function for each time-difference r. The detailed steps are 
given at the end of this subsection. 

The results as well as a comparison with the Gaussian closure method and a direct Monte-Carlo 
simulation are presented in Eigurej^ We can observe that through the proposed approach we are 
able to satisfactorily approximate the correlation function even close to the non-linear regime g = 2 , 
where the Gaussian closure method presents important discrepancies. 
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Figure 5: Correlation functions Cxx and Cxy of the bistable system with system parameters A = 1 , 
ki = —1, and = 1 subjected to Pier son-Moskowitz spectrum, (a) Amplification factor of g = 2 . 
(b) Amplification factor of g = 10. 


Finally, using the computed parameters 7 and closure coefficients, px^x and Px,y^ we can also 
construct the three-dimensional non-Gaussian joint pdf for the response-response-excitation at dif¬ 
ferent time instants. This will be derived based on the three-dimensional Gaussian copula density 
of the following form: 


C(F(x),F(z),G(y)) = 


1 


Vdet R 


exp 


1 

9 

- $-1 (F(*)) ■ 

T 

- $-1 (Fix)) 1 \ 

<p-Uf(z)) 



Z 

[ $-'(G(j/)) J 


[ $-1 (G(2/)) J j 


(32) 


The three-dimensional non-Gaussian joint pdf for the response-response-excitation at different time 
instants can be expressed as follows: 


fx(t),x{t+T),y{t+T){x, z, y) = f(x)f{z)g(y)C (F(x), F{z), G(y)) 


= f{x)f{z)g{y) 


1 


V det R 


exp 



- ^-CFix)) - 

T 

■ - 

\ 

0 

^-CFiz)) 

■ ■ 



z 

\ 

[ ^-CGiy)) \ 


[ <^>-1 (Giy)) \ 



(33) 
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where R represents the 3x3 correlation matrix with all diagonal elements equal to 1: 


R = 


'-'XZ 

^xy C 


'xz ^xy 
1 Czy 

1 


zy 


( 34 ) 


The time dependent parameters Cxz^ Cxy, Czy of the copula function can be found through the 


resolved moments, by expanding the latter as: 

Cxx{r) = jj xzf^(^t),x(t+T),y{t+T){x,z,y)dxdydz = 2J^^Cxz + O , 

Cxy{r) = jJ xyf^(^t^^^(^t+^^^yi^t+^^{x,z,y)dxdydz = 2J^gcxy+o(cly^, 
Cxy(0) = JJ zyf^^t),x(t+T),y{t+T){x,z,y)dxdydz = 2J^gczy + O (^cly^ . 


( 35 ) 

( 36 ) 

( 37 ) 


where, 


J^ = 


j xf{x) erf ^ {2F{x) — 1) dx and ^ ~ J ^9{^) ^^f ^ {2G{x) — 1) dx. 


If necessary higher order terms may be retained in the Taylor expansion although for the present 
problem a linear approximation was sufficient. The computed approximation is presented in Figure 
[^through two dimensional marginals as well as through isosurfaces of the full three-dimensional joint 
pdf. We compare with direct Monte-Carlo simulations and as we are able to observe, the computed 
pdf compares favorably with the expensive Monte-Carlo simulation. The joint statistics using the 
Monte-Carlo approach were computed using 10^ number of samples while the computational cost 
of the MECC method involved the minimization of a three dimensional function. 
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Monte Carlo Simulation 


MECC Method 



Figure 6: Joint pdf fx{t)x{t-\-T)y(t-\-T){^^ computed using direct Monte-Carlo simulation and the 
MECC method. The system parameters are given by A = 1, = — 1, and ks = 1 and the excitation 

is Gaussian following a Pierson-Moskowitz spectrum with g = 10. The pdf is presented through 
two dimensional marginals as well as through isosurfaces. (a)r = 3. (b)r = 10. 


3.2 SDOF bistable oscillator coupled to an electromechanical harvester 

In practical configurations, energy harvesting occurs through a linear electromechanical transducer 


coupled to the nonlinear oscillator 42 60 61 . In this section, we assess how our method performs 
for a bistable nonlinear SDOF oscillator coupled to a linear electromechanical transducer. The 
equations of motion in this case take the form: 


X Xx kix + ksx^ av = y, 

V = Sx, 


(38) 

(39) 


where x is the response displacement, ^ is a stationary stochastic excitation, v is the voltage across 
the load, A is the normalized damping coefficient, ki and ks are the normalized stiffness coefficients, 
a and 6 are the normalized coupling coefficients, and p is the normalized time coefficient for the 
electrical system. All the coefficients except ki are positive. Based on the linearity of the second 
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equation, we express the voltage in an integral form: 


v{t) = 6 f x{C)e ^^d( = 6x{t)^e 

Jo 


( 40 ) 


where * indicates convolution and u{t) represents the Heaviside step function. We then formulate 
the second-order moment equations following a similar approach with the previous section. 




d- 


^2 


-^x{t)y{s) -h X-^x{t)y{s) + kix{t)y{s) -h k 3 x{t)^y{s) + av{t)y{s) = -^y{t)y{s), (41) 


^2 


d 


^2 


-^x{t)x{s) -h X-^x{t)x{s) -h kix{t)x{s) + k 3 x{t)^x{s) + av{t)x{s) = -^y{t)x{s), (42) 


—v{t)v{s) + (3v{t)v{s) = 5—x{t)v{s). (43) 


In this case, we estimate two additional covariance functions, v{t)y{s) and v{t)x{s) before applying 
MECC method: 


Kt)y{s) = S f x{C)y{s)e 
Jo 

= -s)* e~^*'u{t), 


(44) 

(45) 

(46) 

(47) 


where t = t — s is the time difference of two generic time instants t and s. Considering the power 
spectrum, the Fourier transform of the above gives: 


J-{v(t)y{s)} — Syy{Uj) — a . Sxy{(^)- 


jSuj 

P+jLo" 


Similarly, we also obtain for v{t)x{s): 


J^{v{t)x{s)} = Sxx{<^) 


jSu! 


p + juj 




(48) 


(49) 


By applying the previously described closure scheme on equations (41) and (42), we have a linear 
set of moment equations for the second-order statistics: 

■^Cxyir) + X-^Cxyir) + {ki + px,vk3)Cxy{T) + aS-^Cxyir) * e~^^u{t) = -^Cyyir), (50) 
-^Cxx{t) + X-^Cxxir) + {ki + px,xk3)Cxx{T) + aS-^Cxxir) * e~^^u{t) = -^Cxy{-T), (51) 

^C'„„(t) -h PCyxir) = S—Cyx{-T). (52) 
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Using the Wiener-Khinchin theorem, we transform the above equations to the corresponding power 
spectral density equations: 


- c 

+ Px^yks + p ^ ~ 

{(j^) kI Px^xks ~7^^Sxx{p^^ — (j^) ^xyi^)’) 

{-{jUj) + f3}Syy{uj) = - 6{jUj)Syx{(^)- 


(53) 

(54) 

(55) 


These equations allow us to obtain an expression for the power spectral density of the response 
displacement and response voltage in terms of the excitation spectrum: 


Sxx{^) — 


Syvi^^^ - 


{^1 + Px,yk3 — + i(Acj) + ^^}{kl + Px,xk3 “ — j{Xuj) — |z^} 


+ uj‘^}{ki + Px,yk3 - + j{Xuj) + ^^}{kl + Px,xk3 - - j{Xuj) - 


Syy{uj), (56) 


Syy^LU). 


ja6uj ' 


(57) 


Integration of the above equation will give us the variance of the response displacement and voltage: 


/■ 

f 

JQ 


{fcl + P.,yk3 -w^+ jiXuj) + ^}{h + P,,,k3 jiXuj) - 1^} 


0 (jU 


{/32 + co^}{k, + p^^yks -u;^+ j{Xuj) + + p.,.ks - - j{Xoj) - 


Syy{uj)duj, (58) 

Syy{uj)dUJ. 

(59) 


Equation (58) expresses the second order dynamics of the SDOF bistable oscillator coupled with 


an electromechanical harvester, and is the dynamics constraint for this system. We will minimize 
it together with the closure constraints defined by equations (21) and (22): 


■) Px,x-, px,y) — \X^ 


f 


^Syyiiu) 


{kl + px,yk3 — CU^ j{Xcu) + ^^}{kl + px,xk3 — iXp — j{Xoj) — 


duj 


j /x^/(x)erf ^ {2F{x) — 1) dx^ ^ j / j:^/(x) erf ^ {2F{x) — 1) dx^ 

+ "I TWET r + i tETTEEEFW 


/ xf{x) erf ^ {2F{x) — 1) dx 


f xf{x) erf ^ (2F{x) — 1) dx 


(60) 


In Figure we illustrate the variance of the response displacement and the voltage as the 
intensity of the excitation varies for two sets of the system parameters. For both sets of system 
parameters, we observe that for large intensity of the excitation, the MECC method computes the 
response variances (displacement and voltage) very accurately, while the Gaussian closure method 
systematically underestimates them. For lower intensities of the excitation, the response displace¬ 
ment variance computed by the Monte-Carlo simulation presents a non-monotonic behavior with 
respect to q. While the Gaussian closure has very poor performance on capturing this trend, the 
MECC method can still provide a satisfactory approximation of the dynamics. 
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Figure 7: Mean square response displacement and mean square response voltage with respect to 
the amplification factor of Pierson-Moskowitz spectrum for bistable system with two different sets 
of system parameters. Electromechanical harvester parameters are (a = 0.01, /3 = 1, and 6 = 1. (a) 
X = 1, ki = —1, and ks = 1. (b) A = 0.5, ki = —0.5, and ks = 1.0. 

Following similar steps with the previous section, we obtain the covariance functions of the response 
displacement and voltage and the joint pdf in a post-process manner. The results as well as a com¬ 
parison with the Gaussian closure method and the Monte-Carlo simulation are illustrated in Figure 
We can observe that through the proposed approach we are able to satisfactorily approximate 
the correlation function even close to the non-linear regime g = 2, where the Gaussian closure 
method presents important discrepancies. In Figure we illustrate two dimensional marginal pdfs 
as well as isosurfaces of the full three-dimensional joint pdf. We compare with direct Monte-Carlo 
simulations and as we are able to observe, the computed pdf closely approximates the expensive 
Monte-Carlo simulation in statistical regimes which are far from Gaussian. 
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Figure 8: Correlation functions Cxx and of the bistable system with A = 1, = — 1, and /cs = 1 

subjected to Pierson-Moskowitz spectrum. Electromechanical harvester parameters are a = 0.01, 
/3 = 1, and ^ = 1. (a) Amplification factor oi q = 2. (b) Amplification factor of g' = 10. 
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Monte Carlo Simulation 


MECC Method 



Figure 9: Joint pdf fx{t)x{t-\-T)y(t-\-T){^^ computed using direct Monte-Carlo simulation and the 
MECC method. The system parameters are given by A = 1, = — 1, and ks = 1 under Pierson- 

Moskowitz spectrum g = 10. Electromechanical harvester parameters are a = 0.01, P = 1, and 
6 = 1. The pdf is presented through two dimensional marginals as well as through isosurfaces, (a) 
r = 3. (b) r = 10. 

Einally in Eigure[^ we demonstrate how the proposed MECC method can be used to study ro¬ 
bustness over variations of the excitation parameters. In particular, we present the mean square re¬ 
sponse displacement and response voltage estimated for various amplification factors q and frequency- 
varied excitation spectra: 


Sp{uj) = S{uj - cjo), (61) 

where c^o is the perturbation frequency. The comparison with direct Monte-Carlo simulation indi¬ 
cates the effectiveness of the presented method to capture accurately the response characteristics 
over a wide range of input parameters. 
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Monte Carlo Simulation MECC Method Comparison 



Figure 10: Performance comparison (mean square response displacement (a) and voltage (b)) be¬ 
tween Monte-Carlo simulations (100 realizations) and MECC method. Results are shown in terms of 
the amplification factor q and the perturbation frequency cuq of the excitation spectrum (Pierson- 
Moskowitz) for the bistable system with A = 1, = — 1, and /cs = 1. The electromechanical 

harvester parameters are a = 0.01, P = 1, and ^ = 1. 


4 Conclusions 

We have considered the problem of determining the non-Gaussian steady state statistical structure 
of bistable nonlinear vibrational systems subjected to colored noise excitation. We first derived 
moment equations that describe the dynamics governing the two-time statistics. We then combined 
those with a non-Gaussian pdf representation for the joint response-response and joint response- 
excitation statistics. This representation has i) single time statistical structure consistent with 
the analytical solutions of the Fokker-Planck equation, and ii) two-time statistical structure that 
follows from the adoption of a Gaussian copula function. The pdf representation takes the form 
of closure constraints while the moment equations have the form of a dynamics constraint. We 
formulated the two sets of constraints as a low-dimensional minimization problem with respect to 
the unknown parameters of the representation. The minimization of both the dynamics constraint 
and the closure constraints imposes an interplay between these two factors. 

We then applied the presented method to two nonlinear oscillators in the context of vibration 
energy harvesting. One is a single degree of freedom (SDOF) bistable oscillator with linear damp¬ 
ing while the other is a same SDOF bistable oscillator coupled with an electromechanical energy 
harvester. For both applications, it was assumed that the stationary stochastic excitation has a 
power spectral density given by the Pier son-Moskowitz spectrum. We have shown that the pre¬ 
sented method can provide a very good approximation of second order statistics of the system, 
when compared with direct Monte-Garlo simulations, even in essentially nonlinear regimes, where 
Gaussian closure techniques fail completely to capture the dynamics. In addition, we can compute 
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the full (non-Gaussian) probabilistic structure of the solution in a post-process manner. We empha¬ 
size that the computational cost associated with the new method is considerably smaller compared 
with methods that evolve the pdf of the solution since MECC method relies on the minimization 
of a function with a few unknown variables. 

These results indicate that the new method can be a very good candidate when it comes to 
the calculation of the stochastic response for vibrational system with complex potentials as it is 
required in parameter optimization or selection. Future endeavors include the application of the 
presented approach in higher dimensional contexts involving nonlinear energy harvesters and passive 
protection of structures as well as on the development/optimization of structural configurations able 
to operate effectively under intermittent loads |^ . 
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